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I. INTRODUCTION 



Energy transfer (relaxation) phenomena are ubiquitous in nature. At a macroscopic level, 
the phenomenological theory of heat (Fourier law) successfully describes heat transfer and 
energy flow. However, its microscopic origin is still under debate. This is because the phe- 
nomena can contain many-body, multi-scale, nonequilibrium, and even quantum mechanical 
aspects, which present significant challenges to theories addressing energy transfer phenom- 
ena in physics, chemistry and biology jlj]. For example, heat generation and transfer in 
nano-devices is a critical problem in the design of nanotechnology. In molecular physics, it 
i= well known that vibrational energy relaxation (VER) . an e^sentW aspect of any qnanti- 
tative description of chemical reactions |2;]. In the celebrated RRKM theory of an absolute 
reaction rate for isolated molecules, it is assumed that the intramolecular vibrational energy 
relaxation (IVR) is much faster than the reaction itself. Under certain statistical assump- 
tions, the reaction rate can be derived |3| . For chemical reactions in solutions, the transition 
state theory and its extension such as Kramer's theory and the Grote-Hynes theory have 



been developed 

m 



and applied to a variety of chemical systems including biomolecular 



systems [6|. However, one cannot always assume that separation of timescales. It has been 
shown that a conformational transition (or reaction) rate can be modulated by the IVR rate 
. As this brief survey demonstrates, a detailed understanding of IVR or VER is essential 
to study the chemical reaction and conformation change of molecules. 

A relatively well understood class of VER is a single vibrational mode embedded in (vi- 
brational) bath modes. If the coupling between the system and bath modes is weak (or 
assumed to be weak), a Fermi's-golden-rule style formula derived using 2nd order pertur- 
bation theory 8l-ll0| may be used to estimate the VER rate. However, the application of 
such theories to real molecular systems poses several (technical) challenges, including how 
to choose force fields, how to separate quantum and classical degrees of freedom, or how to 
treat the separation of timescales between system and bath modes. Multiple solutions have 
been proposed to meet those challenges leading to a variety of theoretical approaches to the 



treatment of VER 



11 



Iq . These works using Fermi's golden rule are based on quantum 



mechanics and suitable for the description of high frequency modes (more than thermal 
energy ^ 200 cm~^), on which nonlinear spectroscopy has recently focused 17H20|. 

In this chapter, we summarize our recent work on VER of high frequency modes in 
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biomolecular systems. In our previous work, we have concentrated on the VER rate and 



mechanisms for proteins 



2l| . Here we shall focus on the time course of the VER dynamics. 



We extend our previous Markovian theory of VER to a non-Markovian theory applicable 



to a broader range of chemical systems 



22 



23|. Recent time- resolved spectroscopy can 



detect the time course of VER dynamics (with femto-second resolution), which may not be 
accurately described by a single time scale. We derive new formulas for VER dynamics, and 
apply them to several interesting cases, where comparison to experimental data is available. 

This chapter is organized as follows: In Sec. Ull we briefly summarize the normal mode 
concepts in protein dynamics simulations, on which we build our non-Markovian VER theory. 
In Sec. mil we derive VER formulas under several assumptions, and discuss the limitations 
of our formulas. In Sec. IIVI we apply the VER formulas to several situations: the amide 
I modes in solvated A^-methylacetamide and cytochrome c, and two in-plane modes (z/4 
and z/7 modes) in a porphyrin ligated to imidazol. We employ a number of approximations 



in describing the potentia 
the empirical CHARMM 



energy surface on which the dynamics takes place, including 



24( 1 force field and density functional calculations 25| for the 



small parts of the system (A^-methylacetamide and porphyrin). We compare our theoretical 
results with experiment when available, and find good agreement. We can deduce the VER 
mechanism based on our theory for each case. In Sec. |Vl we summarize and discuss the 
further aspects of VER in biomolecules and in nanotechnology (molecular devices). 



II. NORMAL MODE CONCEPTS APPLIED TO PROTEIN DYNAMICS 



Normal mode provides a powerful tool in exploring molecular vibrational dynamics 26 1 
and may be applied to biomolecules as well [271. The first normal mode calculations for 
a protein were performed for BPTI protein |28| . Most biomolecular simulation softwares 



support the calculation of normal modes 



However, the calculation of a mass- 



weighted hessian Kij^ which requires the second derivatives of the potential energy surface, 
with elements defined as 



1 



9V 



(1) 



can be computational demanding. Here rrii is the mass, Xj is the coordinate, and V is the 
potential energy of the system. Efficient methods have been devised including torsional angle 
normal mode 311], block normal mode 32|, and the iterative mixed-basis diagonalization 



(DIMB) methods 33|] etc. An alternative direction for the efficient calculations of a hessian 



34| or Gaussian network 



models. From 



is to use coarse-grained models such as elastic 
normal mode analysis (or instantaneous normal mode analysis [36|), the frequencies, density 
of states, and normal mode vectors can be calculated. In particular, the latter quantity is 
important because it is known that the lowest eigenvectors may describe the functionally 
important motions such as large-scale conformational change, a subject that is the focus of 
another chapter in this volume js^]. 

There is no doubt as to the usefulness of normal mode concepts. However, for molecular 
systems, it is always an approximate model as higher-order nonlinear coupling and intrinsic 
anharmonicity become essential. To describe energy transfer (or relaxation) phenomena 
in a protein, Moritsugu, Miyashita, and Kidera (MMK) introduced a reduced model using 
normal modes with 3rd and 4th order anharmonicity [38], C'^j^ and C^f^,^, respectively. 



v({*}) = Et''' + 5!^^' 

k klm 



(3) , 
klmlkQlQm + 



klmn 



(4) 

klmnlklllmQn 



with 



C, 



(3) 

klm 



C, 



(4) 

klmn 



dqkdqidqm' 



(2) 

(3) 
(4) 



dqkdqidqmdqn 

where qk denotes the normal mode calculated by the hessian Kij and Uk is the normal mode 



frequency. Classical (and harmonic) Fermi resonance 39| is a key ingredient in the MMK 
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FIG. 1: Left: The excited eigenvector depicted by arrows in myoglobin. Right: Classical simulation 
of mode specific energy transfer in myoglobin at zero temperature. [Reprinted with permission from 
K. Moritsugu, O. Miyashita and A. Kidera, Phys. Rev. Lett. 85, 3970 (2000). Copyright ©2009 
by the American Physical Society.] 
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theory of energy transfer derived from observations of all-atom simulations of myoglobin at 
zero temperature (see Fig. [1]). 

At finite temperature, non-resonant effects become important and clear interpretation 
of the numerical results becomes difficult within the classical approximation. Nagaoka and 
coworkers 40| identified essential vibrational modes in vacuum simulations of myoglobin 
and connected these modes to the mechanism of "heme cooling" explored experimentally by 
Mizutani and Kitagawa [isl. Contemporaneously, nonequilibrium MD simulations of sol- 
vated myoglobin carried out by Sagnella and Straub provided the first detailed and accurate 
simulation of heme cooling dynamics [41]. That work provided support for the conjecture 
that the motion similar to those modes identified by Nagaoka play an important role in 
energy flow pathways. 

Nguyen and Stock explored the vibrational dynamics of the small molecule, N- 
methylacetamide (NMA), often used as a model of the peptide backbone [42|. Using 
nonequilibrium MD simulations of NMA in heavy water, VER was observed to occur on 
a pico-second time scale for the arnide I vibrational mode (see Fig. |2]). They used the in- 
stantaneous normal mode concept ^] to interpret their result and noted the essential role 
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FIG. 2: Nonequilibrium MD simulation of energy flow from the excited amide I mode in A'^- 
methylacetamide in heavy water. See also Fig. [3l [Reprinted with permission from P.H. Nguyen 
and G. Stock, J. Chem. Phys. 119, 11350 (2003). Copyright ©2009 by the American Institute of 
Physics.] 
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of anharmonic coupling. Leitner also used the normal mode concept to describe energy 
diffusion in a protein, and found an interesting link between the anomalous heat diffusion 



and the geometrical properties of a protein 



43|. 



In terms of vibrational spectroscopy, Gerber et al. calculated the anharmonic frequencies 



in BPTI, within the VSCF level of theory 4J], using the reduced model ([2]). Yagi, Hirata 



and Hirao refined this type of anharmoic frequency calculation for large molecular systems 



with more efficient methods 



451], appropriate for applications to biomolecules such as DNA 
base pair [46| . Based on the reduced model ([2]) with higher-order nonlinear coupling, Leitner 
also studied quantum mechanical aspects of VER in proteins, by employing the Maradudin- 
Fein theory based on Fermi's "Golden Rule" [ij]. Using the same model, Fujisaki, Zhang, 
and Straub focused on more detailed aspects of VER in biomolecular systems, and calculated 
the VER rate, mechanisms or pathways, using their non-Markovian perturbative formulas 
(described below). 

As this brief survey demonstrates, the normal mode concept is a powerful tool that 
provides significant insight into mode specific vibrational dynamics and energy transfer in 
proteins. 

III. DERIVATION OF NON-MARKOVIAN VER FORMULAS 

We have derived a VER formula for the simplest situation, a one- dimensional relaxing 



oscillator coupled to a "static" bath j22|]. Here we extend this treatment to two more general 
directions: (a) multi-dimensional relaxing modes coupled to a "static" bath and (b) a one- 
dimensional relaxing mode coupled to a "fluctuating" bath 

A. Multi-dimensional relaxing mode coupled to a static bath 

We take the following time-independent Hamiltonian 

n = H% + Hb + V^ (5) 

= ul + {v)B + nB + v'^ -{v)b (6) 
= ns + nB + v (7) 
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where 



Hs = n% + (V)b, (8) 
V = V^- {V)b. (9) 



In previous work [22], we have considered only a single one- dimensional oscillator as the 
system. Here we extend that treatment to the case of an A^5-dimensional oscillator system. 
That is, 

i=l ^ ^ 

= E (f + T«°) • 

Ns 

V = (12) 

1=1 

where V{{qi}) is the interaction potential function between Ns system modes which can be 
described by, for example, the reduced model, Eq. ([2]). The simplest case V{{qi}) = is triv- 
ial as each system mode may be treated seperately within the perturbation approximation 
for V. 

We assume that \k) is & certain state in the Hilbert space spanned by l-Ls- Then the 
reduced density matrix is 

{psUnit) = (m|e-^«^*/^Tr5{p(t)}e^«-*/^|n) (13) 

where tilde denotes the interaction picture. Substituting the time-dependent perturbation 
expansion 

1 /■* ~ 



m = p(o) + ^ J dt'[vit'),pio)] 

ft rt' 



into the above, we find 



ly^ldt'l^ rft"[V(t'),[V(t"),p(0)]] + --- (14) 



{Ps)mn{t) ^ (psYlit) + {Ps)l^lit) + ips)ilit) + ■ ■ ■ , (15) 
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where 



{m{-t)\psmn{-t)) = {m\ps{t)\n) (16) 

t rt' 



iPsflit) = dt' rft"(m|e-^«-*/^TrB{[V(t'),[V(t"),p(0)]]}e^««*/^|n) 

1 



W 11 fo ^^""^"^^ ' W)qAt")Ps{Q) 
-qAt'')psmrmn{-t)){5Ht')5^At''))B 

+ 7w f^^' f dt"Y,{^i-mpsmAt"Mt') 

i'!''!') Jo Jo 

-g.(t')P5(0)g,(t")]|n(-t))(5J-,(n5J-.(t'))i^. (17) 

Here we have defined \m{t)) = e~*^'^*''^|m) and taken {ps)rnn{t) = 0. Recognizing that we 
must evaluate expressions of the form 

Rran;ij{t;t',t") = {m{-t)\[qi{t%{t") psiO)\n{-t)) , 

-{m{-t)\q,{npsm{t')]\n{-t)) (18) 
a,it',t") = {5F,{t')5F,{t"))B (19) 

and their complex conjugates, R^,^.ij{t]t',t"),C*j{t',t"), the 2nd order contribution can be 
written 

+R:^.^,^it-,t',t'')c:^it',t'')]. (20) 

We can seperately treat the two terms. Assuming that we can solve 7^5 1 a) = Ea\a), we find 

{t;t',t") = y^Xm\a){qi)ab{qj)bc{ps)cd{d\n) 

abed 

y^^-i{Ea-Ea)t-i(Et-Ea)t'-i(E,-Et)t" 

- '^{m\a){qj)abips)bc{qi)cd{d\n) 

abed 

^ ^-i{Ea-Ea)t-i{Ea-Ec)t'-i(Eh-Ea)t" ^21) 

For the bath-averaged term, we assume the following force due to third-order nonlinear 



coupling of system mode i to the normal modes, a and /3, of the bath 21 1 



SJ^iiiqa}) = Cia^iqaq^ " (qaq^)) (22) 

a,/3 
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and we have 21 1 



a,{t', t") = BTf^', t") + Bi^t', t") + t") 



(23) 



with 



R±-{t',t") 



a, 13 



i{LUa+Ulfl){t'-t") 



a, 13 



where 



D. 



al3;ij 



Ciaj3Cja/3 

and Ha is the thermal population of the bath mode a. 



(24) 

(25) 
(26) 

(27) 



This formula reduces to our previous result for a one-dimensional system oscillator 22 1 
when Ns = 1 and all indices are suppressed. Importantly, this formula can be ap- 

plied to situations where it is difficult to define a "good" normal mode to serve as a one- 
dimensional relaxing mode, as in the case of the CH stretching modes of a methyl group 21|. 
However, expanding to an A^^ dimensional system adds the burden of solving the multidi- 
mensional Schrodinger equation 7^5 1 a) = Ea\a). To address this challenge we may employ 
vibrational self-consistent field (VSCF) theory and its extensions developed by Bowman 
and coworkers 48|] implemented in MULTIMODE program of Carter and Bowman 49j or 
in the SINDO program 50(] of Yagi and coworkers. As in the case of our previous theory 
of a one- dimensional system mode, we must calculate A^5-tiple 3rd order coupling constants 
Ciapii = 1, 2, ...Ns) for all bath modes a and (3. 



B. One-dimensional relaxing mode coupled to a fluctuating bath 

We start from the following time- dependent Hamiltonian 

-Hit) = nl{t) + V,B{t) + V\t) (28) 
= nl{t) + (V(t))B + UBit) + V\t) - {V{t))B (29) 

= Us{t) + nB{t) + v{t) (30) 
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where 

Us{t) ^ U%{t) + {V{t))B, 
V{t) ^ V\t) - {V{t))B 

with the goal of solving the time-dependent Schrodinger equation 

^^^^ = [ns{t) + nB{t) + v{tm{t)) = [no{t) + v{tm{t))- 

By introducing a unitary operator C/o(i) = Us{t)UB{t) 

i^jU^it) = -Homoit), 

ihfusit) = nsimsit), 

ih^UBit) = nB{t)UB{t), 
at 

we can derive an "interaction picture" von Neumann equation 
where 

v{t) = uiit)vmo{t), 
m = ui{t)p{t)Uo{t). 

We assume the simple form of a harmonic system and bath, but allow fluctuations 
system and bath modes modeled by time-dependent frequencies 

Hsit) = noos{t){alas + 1/2), 
V-Bit) = ^huJo.{t){aiao, + 1/2). 

a 

The unitary operators generated by these Hamiltonians are 

and the time evolution of the annihilation operators is given by 

Ul{t)asUs{t) = ase-^^odro^sir)^ 

10 



To simplify the evaluation of the force autocorrelation function, we assume that the temper- 
ature is low or the system mode frequency is high as a justification for the approximation. 
Substituting the above result into the force autocorrelation function calculated by the force 
operator, Eq. fl22]) . we find 

(ST(t')SJ'(t")) ~ Csal3it')Csal3it") 

y^(,-ilBcp{t')-ec(l{t")] ^4g^ 

where 

Qsit) = [ drusir), (47) 

Qa^it) = I dT[tUo^{T)+tUp{T)]. (48) 

Jo 

Substituting this approximation into the perturbation expansion Eqs. ( TT^ . ( ITB]) . and (IT7|1 . 
we obtain our final result 

C Sal3{'t')C Sapit") 



^us{t')ujc,{t')up{t')ujs{t'')ujc,{t'')ujp{V') 
X COS {e5(t') - Qapit') - Qs{t") + Qapit")} (49) 



which provides a dynamic correction to the previous formula [22j. The time- dependent 
parameters uJs{t),uJa{t),Csai3{f) may be computed from a running trajectory using instan- 
taneous normal mode analysis jsil. This result was first derived by Fujisaki and Stock 4^], 
and applied to the VER dynamics of A^-methylacetamide as described below. This correction 
eliminates the assumption that the bath frequencies are static on the VER timescale. 

For the case of a static bath, the frequency and coupling parameters are time-independent 
and this formula reduces to the previous one-dimensional formula (when the off-resonant 
terms are neglected) 



22 



/ N / ,N h ^ Cg^p 1 - COS (^5 - - /Kn^ 

Note that Bakker derived a similar fiuctuating Landau- Teller formula in a different man- 



ner 



5l| . It was successfully applied to molecular systems by Sibert and coworkers [52 1. 
However, the above formula differs from Bakker's as (a) we use the instantaneous normal 
mode analysis to parameterize our expression, and (b) we do not take the Markov limit. 
Our formula can describe the time course of the density matrix as well as the VER rate. 
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One further point is that we use the cumulant-type approximation to calculate the dy- 
namics. When we calculate an excited state probability, we use 



(P5)ii(t) = 1 - (psMt) ~ exp{-(p5)oo(t)}. (51) 

Of course, this is vahd for the initial process {{ps)oo(t) ^ l); but, at longer time scales, we 
take {ps)ii{t) ~ exT>i — (ps)nn(t)} because the naive formula {ps)ii(t) = 1 — (ps')oo(^) can be 



negative, which is unphysical 47 1. 



C. Limitations of the VER formulas and comments 

There are several limitations to the VER formulas derived above. The most obvious is 
that they are 2nd order perturbative formulas and rely on a short-time approximation. As 
far as we know, however, there exists no non-perturbative quantum mechanical treatment 
of VER applicable to large molecular systems. It is prohibitive to treat the full molecular 



dynamics quantum mechanically 



53 1 for large molecules. Moreover, while there exist several 



mixed quantum-classical methods 11| that may be applied to the study of VER, but there is 



no guarantee that such approximate methods work better than the perturbative treatment 



54|. 



Another important limitation is the adaptation of a normal mode basis set, a natural 
choice for molecular vibrations. Because of the normal mode analysis, the computation 

there is a 



can be burdensome. When we employ instantaneous normal mode analysis |; 
concern about the imaginary frequency modes. For the study of high frequency modes, this 
may not be significant. However, for the study of low frequency modes, the divergence of 
quantum (or classical) dynamics due to the presence of such imaginary frequency modes is 
a significant concern. For the study of low frequency modes, it is more satisfactory to use 



other methods that do not rely on normal mode analysis such as semiclassical methods 55 1 



or path integral methods 



56] 



We often use "empirical" force fields to describe quantum dynamics. However, it is well 
known that the force fields underestimate anharmonicity of molecular vibrations [57]. It 
is often desirable to use ab initio potential energy surfaces. However, that more rigorous 
approach can be demanding. Lower levels of theory can fail to match the accuracy of some 
empirical potentials. As a compromise, approximate potentials of intermediate accuracy. 
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such as QM/MM potentials [58|, may be appropriate. We discuss this question further in 
Sec. WJ\ and in Sec. HVCl 



IV. APPLICATIONS OF THE VER FORMULAS TO VIBRATIONAL MODES IN 
BIOMOLECULES 

We study the quantum dynamics of high frequency modes in biomolecular systems using 
a variety of VER formulas described in Sec. IIIII The application of a variety of theoreti- 
cal approaches to essential VER processes will allow for a relative comparison of theories 
as well as the absolute assessment of theoretical predictions compared with experimental 
observations. In doing so, we address a number of fundamental questions. What are the 
limitations of the static bath approximation for fast VER in biomolecular systems? Can the 
relaxation dynamics of a relaxing amide I vibration in a protein be accurately modeled as 
a one-dimensional system mode coupled to a harmonic bath? Can the "fluctuating bath" 
model accurately capture the system dynamics when the static picture of normal modes is 
not "good" on the timescale of the VER process? In Sec. IIV Al and lIVBt our main focus 
is the VER of excited amide I modes in peptides or proteins. In Sec. IIV Ct we study some 
vibrational modes in porphyrin ligated to imizadol, which is a mimic of a heme molecule in 
heme-proteins including myoglobin and hemoglobin. 



A. iV-methylacetamide 



iV-methylacetamide (NMA) is a well-studied small molecule (CH3-CO-NH-CH3) that 
serves as a convenient model of a peptide bond structure (-CO-NH-) in theory and experi- 
ment. As in other amino acids, there is an amide I mode, localized on the CO bond stretch, 
which is a useful "reporter" of peptide structure and dynamics when probed by infrared 
spectrocopy. Many theoretical and experimental studies on amide I and other vibrational 
modes (amide II, III) have characterized how the mode frequencies depend on the local sec- 



ondary structure of peptides or proteins [59 



and polarizability of these modes, see Refs. 



60 



15 



16 



or t 



61 



le accurate description of frequencies 



65| . The main focus of these works is 



the frequency sensitivity of amide modes on the molecular configuration and environment. 
In this case, the amide mode frequencies are treated in a quantum mechanical way, but the 
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configuration is treated classically. With a focus on interpreting mode frequency shifts due 
to configuration and environment, mode-coupling between amide modes and other modes is 
often neglected. As we are mainly interested in VER or IVR dynamics of these modes, an 
accurate treatment of the mode coupling is essential. 

Recent theoretical development of IVR dynamics in small molecules is summarized in 53 1. 
Leitner and Wolynes {?!] utilized the concept of local random matrix to clarify the quantum 
aspects of such dynamics. The usefulness and applications of their approach are summarized 
in 12(1 as well as in this volume 13|]. However, these studies are focused on isolated molecules, 



whereas our main interest is in exploring quantum dynamics in a condensed phase. We take 
a step-by-step hierarchical approach. Starting from the isolated NMA molecule, we add 
several water molecules to form NMA-water clusters, and finally treat the condensed phase 
NMA- water system (see Fig [3]). With increasing complexity of our dynamical model, the 
accuracy of our theory, including the quality of the potential energy surface and the accuracy 
of the quantum dynamics must diminish. As such, a pricipal focus of our account is a careful 
examination and validation of our procedures through comparison with accurate methods 
or experiment. 




FIG. 3: Representation of three models employed for the study of VER dynamics in A^- 
methylacetamide (NMA). (a) NMA, (b) NMA with three solvating water, (c) NMA with first 
solvation shell derived from simulations in bulk water. [Reprinted with permission for (a) and 
(b) from Y. Zhang, H. Fujisaki, and J.E. Straub, J. Phys. Chem. A 113, 3051 (2009). Copyright 
©2009 by the American Chemical Society.] [Reprinted with permission for (c) from H. Fujisaki 
and G. Stock, J. Chem. Phys. 129,134110 (2008). Copyright ©2009 by the American Institute of 
Physics.] 
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1. N -methylacetamide in vacuum 



In our studies of isolated NMA 



66 



671], we have employed both accurate potential en- 



ergy surface and accurate quantum dynamics methods to explore the timescale and mecha- 
nism of VER. From the anharmonic frequency calculations and comparison to experiment 
we concluded that B3LYP/6-31G(d) is a method of choice for computation of the 



electronic ground state potential surface, considering both accuracy and feasibility. For 



other 
Refs. 



;reatments at differing levels of theory of quantum chemical calculation on NMA, see 



57|, 



69 



70|. After the construction of an accurate potential surface, there are several 



tractable approaches for treating the quantum dynamics for this system. The most accurate 
is the vibrational configuration interaction (VCI) method based on vibrational self-consistent 



field (VSCF) basis sets (see Refs. j48|, |49|, |66|, |67| for details). We employed the Sindo code 
developed by Yagi 50| . The numerical results for the VCI calculation are shown in Fig. HJ 
and compared to the prediction based on the perturbative formula Eq. fl50l) and classical 
calculations as done in j42j]. Both approximate methods seem to work well, but there are 
caveats. The perturbative formula only works at short time scales. There is ambiguity for 
the classical simulation regarding how the zero point energy correction should be included 
(see Stock's papers |71|). The main results for a singly deuterated NMA (NMA-rfi) are (1) 
the relaxation time appears to be sub ps, (2) as NMA is a small molecule, there is a recurrent 
phenomenon, (3) the dominant relaxation pathway involves three bath modes as shown in 
Fig. [5l and (4) the dominant pathways can be identified and characterized by the following 



Fermi resonance parameter 



66 



67| 



mv\f) 



AE 



oc 



C, 



Ski 



h 



h 



2ujq V 2a;, 



h 

2^1 



(52) 



where {i\AV\f) is the matrix element for the anharmonic coupling interaction, and AE = 
h{us — — u^i) is the resonance condition (frequency matching) for the system and two 
bath modes. Both the resonant condition (AE) and anharmonic coupling elements (Cski) 
play a role, but we found that the former affects the result more significantly. This indi- 
ates that, for the description of VER phenomena in molecules, accurate calculation of the 
harmonic frequencies is more important than the accurate calculation of anharmonic cou- 
pling elements. This observation is the basis for the development and application of the 
multi- resolution methods for anharmonic frequency calculations [3]. 
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2. N-methylacetamide/Water cluster 



We next examine a somewhat larger system, NMA in a water cluster 72|, an interesting 
and important model system for exploring the response of amide vibrational modes to "sol- 
vation" 62|]. The system size allows for an ab initio quantum mechanical treatment of the 



potential surface at a higher level of theory, B3LYP/aug-cc-pvdz, relative to the commonly 
employed B3LYP/6-31G(d). The enhancement in level of theory significantly improves the 
quality of the NMA-water interaction, specifically the structure and energetics of hydrogen 
bonding. Since there are at most three hydrogen bonding sites in NMA, it is natural to 
configure three water molecules around NMA as a minimal model of "full solvation." NMA- 
water hydrogen bonding causes the frequency of the amide I mode to red-shift. As a result, 
the anharmonic coupling between the relaxing mode and the other bath modes will change 
relative to the case of the isolated NMA. Nevertheless, we observe that the VER time scale 
remains sub ps as is the case for isolated NMA (Fig. [6]). Though there are intermolecular 
(NMA-water) contributions to VER, they do not significantly alter the VER timescale. An- 
other important finding is that the energy pathway from the amide I to amide II mode is 
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FIG. 4: Left: Time evolution of the energy content of the initially excited amide I mode as well 
as all the remaining modes of A^-methylacetamide. Quantum (red lines) and classical (green lines) 
calculations obtained at the DFT/B3LYP level of theory are compared. Right: Comparison of 
the VCI calculation (red) with the result of the perturbative calculation (green) for the reduced 
density matrix. [Reprinted with permission from H. Fujisaki, K. Yagi, J.E. Straub, and G. Stock, 
Int. J. Quant. Chem. 109, 2047 (2009). Copyright ©2009 by Wiley InterScience.] 
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"open" for the NMA-water cluster system. This result is in agreement with experimental re- 



sults by Tokmakoff and coworkers 



74j ] and recent theoretical investigation 16| . Comparison 



between singly (NMA-rfi) and fully (NMA-rfy) deuterated cases shows the VER time scale 
become somewhat longer for the case of NMA-rfy (Fig. [H]). We also discuss this phenomenon 
below in the context of the NMA/solvent water system. 



3. N -methylacetamide in water solvent 



58|. We 



Finally we consider the condensed phase system of NMA in bulk water [22 
attempt to include the full dynamic effect of the system by generating many configurations 
from molecular dynamics simulations and using them to ensemble-average the results. Note 
that in the previous examples of isolated NMA and NMA/water clusters, only one configu- 
ration at a local minimum of the zero temperature ground state potential surface was used. 
On the other hand, the potential energy function used is not so accurate as in the previous 
examples as it is not feasible to include many water molecules at a high level of theory. 
We have used the CHARMM force field to calculate the potential energy and to carry out 
molecular dynamics simulations. 



All simulations were performed using the CHARMM simulation program package 24 1 
and the CHARMM22 all-atom force field ^] was employed to model the solute NMA-c?! 



and the TIP3P water model 



76| with doubled hydrogen masses to model the solvent D2O. 



We also performed simulations for fully deuterated NMA-c/y. The peptide was placed in 






FIG. 5: The dominant bath vibrational modes coupled to the amide I mode calculated on the 
B3LYP/6-31G(d) potential energy surface. [Reprinted with permission from H. Fujisaki, K. Yagi, 
K. Hirao, and J.E. Straub, Chem. Phys. Lett. 443, 6 (2007). Copyright ©2009 by Elsevier.] 
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a periodic cubic box of (25.5 A)^ containing 551 D2O molecules. All bonds containing 
hydrogens were constrained using the SHAKE algorithm {t^]. We used a 10 A cutoff with a 
switching function for the nonbonded interaction calculations. After a standard equilibration 
protocol, we ran a 100 ps NVT trajectory at 300 K, from which 100 statistically independent 
configurations were collected. 

We applied the simplest VER formula Eq. (150|1 |22| as shown in Fig. [3 We truncated 
the system including only NMA and several water molecules around NMA with a cutoff 
distance, taken to be 10 A. For reasons of computational feasibility, we only calculated the 
normal modes and anharmonic coupling elements within this subsystem. 

A number of important conclusions were drawn from these calculations. 
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(1) The inclusion of "many" solvating water molecules induces the irreversible decay of 
the excess energy as well as the density matrix elements (population). The important 
observation is that the VER behavior does not severely depend on the cutoff distance (if 
it is large enough) and the cutoff frequency. The implication is that if we are interested 
in a localized mode such as the amide I mode in NMA, it is enough to use a NMA/water 
cluster system to totally describe the initial process of VER. In a subsequent study, Fujisaki 
and Stock used orily 16 water molecules surrounding NMA (hydrated water) and found 



reasonable results 



47|. 



(2) Comparison of the two isolated NMA calculations suggests that the CHARMM force 
field works well compared with results based on DFT calculations. This suggests that the 
use of the empirical force field in exploring VER of the amide I mode may be justified. 

(3) There is a classical limit of this calculation j22j] , which predicts a slower VER rate close to 
Nguyen-Stock's quasi-classical calculation j^]. This finding was explored further by Stock 
781], who derived a novel quantum correction factor based on the reduced model, Eq. ([2]). 

In these calculations, many solvating water configurations were generated using MD sim- 
ulations. As such, information characterizing dynamic fluctuation in the environment is 

n 

ignored. Fujisaki and Stock further improved the methodology to calculate VER [47] by 
taking into account the dynamic effects of the environment through the incorporation of 
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time-dependent parameters, such as the normal mode frequencies and anharmonic couphng, 
derived from the MD simulations as shown in Fig. |8l Their method is described in Sec. IIII B[ 
and was applied to the same NMA/solvent water system. As we are principally concerned 
with high frequency modes, and the instantaneous normal mode frequencies can become 
unphysical, we adopted a partial optimization strategy. We optimized the NMA under the 
influence of the solvent water at a fixed position. (For a different strategy, see Ref. {tqI.) 
The right side of Fig. |8] shows the numerical result of the optimization procedure. Through 
partial optimization, the fluctuations of the parameters become milder than the previous 
calculations that employed instantaneous normal modes. 

The results of numerical calculations based on Eq. (1491) are shown in Fig. O We see that 
both partial optimization and dynamical averaging affect the result. The "dynamic" formula, 
Eq. fH9|) . leads to smaller fluctuations in the results for the density matrix. Apparently, 
dynamic averaging smoothens the resonant effect, stemming from the frequency difference 
in the denominator of Eq. ( l50|) . For the NMA/solvent water system, the time-averaged value 
of the Fermi resonance parameter, Eg. ( |52l) . can be utilized to clarify the VER pathways 



as in the case of isolated NMA [47l|. It was shown that the hydrated water (the number 
of waters is 16) is enough to fully describe the VER process at the initial stage (^ 0.5 
ps). The predictions of the VER rates for the two deuterated cases, NMA-rfi and NMA-dj, 
are in good agreement with experiment and also with the NMA/Water cluster calculations 
[3]. Though the dynamic effect is modest in the case of the NMA/solvent water system, 
the dynamic formula is recommended when variations of the system parameters due to the 
fluctuating environment must be taken into account. 



B. Cytochrome c in water 



The pro tein, cytochrome c, has been used in experimental and theoretical studies of VER 
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SOMSSl]. Importantly, spectroscopy and simulation have been used to explore the time 



scales and mechanism of VER of CH stretching modes 
amide I modes in cytochrome c 



21 



. Here we examine VER of 



2l| which (a) employed 



Distinct from previous studies 
a static local minimum of the system, we use the dynamical trajectory; (b) in the previous 
study, the water degrees of freedom were excluded, whereas in this study some hydrating 
water has been taken into account. 
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We used the trajectory of cytochrome c in water generated by Bu and Straub |85|. To 
study the local nature of the amide I modes and the correspondence with experiment, we 
isotopically labeled four specific CO bonds, typically C^^O^^ as C^^O^^. In evaluating the 
potential energy in our instantaneous normal mode analysis, we truncated the system with 
an amide I mode at the center using a cutoff (~ 10 A), including both protein and water. 
Following INM analysis, we used Eq. (150|) to calculate the time course of the density matrix. 
The predicted VER is single exponential in character with time scales that are subpicosecond 
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FIG. 8: Time evolution of the vibrational dynamics of NMA in D2O obtained from instantaneous 
normal mode analysis with (right) and without (left) partial energy minimization. Shown are 
(upper panels) the frequency mismatch Ausafsit) = <^sit) — ^a{t) — ^pit), for several resonant 
bath mode combinations, and (lower panels) the corresponding third-order anharmonic couplings, 
Csapit)- [Reprinted with permission from H. Fujisaki and G. Stock, J. Chem. Phys. 129,134110 
(2008). Copyright (c)2009 by the American Institute of Physics.] 
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with relatively small variations induced by the different environments of the amide I modes 
(see Fig. [10] and Table I in 23| for numerical values of the VER timescales). To identify 
the principal contributions to the dependence on the environment we examined the VER 
pathways and the roles played by protein and water degrees of freedom in VER. Our first 
conclusion is that, for the amide I modes buried in the protein (a-helical regions), the 
water contribution is less than for the amide I modes exposed to water (loop regions). This 
finding is important because only a total VER timescale is accessible in experiment. With 
our method, the energy flow pathways into protein or water can be clarified. 

Focusing on the resonant bath modes, we analyzed the anisotropy of the energy flow, 
as shown in Fig. [TTl where the relative positions of bath modes participating in VER are 
projected on the spherical polar coordinates {6, (p) centered on the CO bond involved in the 
amide I mode, which represents the principal z-axis (see Fig. 1 in 23|). The angle dependence 
of the energy flow from the amide I mode to water is calculated from the normal mode 
amplitude average, and not directly related to experimental observables. As expected, energy 
flow is observed in the direction of solvating water. However, that distribution is not spatially 
isotropic and indicates preferential directed energy flow. These calculations demonstrate the 
power of our theoretical analysis in elucidating pathways for spatially directed energy flow 
of fundamental importance to studies of energy flow and signaling in biomolecules and the 
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FIG. 9: VER calculations of amide I mode population P{t) of NMA with use of instantaneous 
normal mode analysis and partial energy minimization. Shown are results from (left) the inho- 
mogeneous averaging approximation and (right) dynamical averaging. [Reprinted with permission 
from H. Fujisaki and G. Stock, J. Chem. Phys. 129,134110 (2008). Copyright ©2009 by the 
American Institute of Physics.] 
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optimal design of nanodevices (see summary and discussion for more detail). 



C. Porphyrin 



Our final example is a modified porphyrin [86:1]. We have carried out systematic studies 
of VER in the porphyrin-imidazole complex, a system that mimics the active site of the 
leme protein myoglobin (Mb). The structure of myoglobin was first determined in 1958 



871]. Experimental and computational studies exploring the dynamics of myoglobin led to 



the first detailed picture of how fluctuations in a protein structure between a multitude 
of "conformational substates" supports protein function [88]. Time resolved spectroscopic 
studies [l^ coupled with computational studies have provided a detailed picture of time 
scale and mechanism for energy flow in myoglobin and its relation to function. Karplus and 
coworkers developed the CHARMM force field 2J] for heme and for amino acids for the 
study of myoglobin, and a particular focus on the dissociation and rebinding of ligands such 





FIG. 10: Left: 81st and 84th residues of cytochrome c in a loop region. Right: 93rd and 97th 
residues of cytochrome c in an a hehcal region. The cartoon represents the protein using a licorice 
model to identify the four residues. (The water molecules are excluded for simplicity.) [Reprinted 
with permission from H. Fujisaki and J.E. Straub, J. Phys. Chem. B 111, 12017 (2007). Copyright 
(c)2009 by the American Chemical Society.] 
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as CO, NO, and O2 90|]. The empirical force field appears to provide an accurate model 
of heme structure and fluctuations. However, we have less confidence in the accuracy of 
anharmonicity essential to detailed mode coupling on the force field and the dependence 
on spin state, which is important to the proper identification of the electronic ground state 
potential energy surface. 

We carried out ab initio calculations for a heme-mimicking molecule, iron-porphine ligated 
to imidazole, abbreviated as FeP-Im. See Fig. [12] for the optimized structure. We employed 
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FIG. 11: Angular excitation functions for the resonant normal modes of water for the (a) 81st, (b) 
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the UB3LYP/6-31G(d) level of theory as in the case of the isolated NMA |66|, |67|, but 
carefully investigated the spin configurations. We identified the quintuplet (5* = 2) as the 
electronic ground state, in accord with experiment. Our study of VER dynamics on this 
quintuplet ground state potential energy surface (PES) is summarized here. Additional 
investigations of the VER dynamics on the PES corresponding to other spin configurations 
as well as different heme models are described elsewhere |86|. 

A series of elegant pioneering experimental studies have provided a detailed picture of 
the dynamics of the 1/4 and z/7 modes, in-plane modes of the heme (see Fig. [T3|l . following 
ligand photodissociation in myoglobin. Using time-resolved resonance Raman spectroscopy, 



Mizutani and Kitagawa observed mode specific excitation and relaxation [18|, |89|]. Interest- 
ingly, these modes decay on different time scales. The VER time scales are ~ 1.0 ps for 
the z/4 mode and ~ 2.0 ps for the z/7 mode. Using a sub-lO-fs pulse. Miller and co-workers 
extended the range of the coherence spectroscopy up to 3000 cm~^ '9l|. The heme mode 
was found to be most strongly excited following Q band excitation. By comparing to the 
deoxyMb spectrum, they demonstrated that the signal was derived from the structural tran- 
sition from the six-coordinate to the five-coordinate heme. Less prominent excitation of the 
1^4 mode was also observed. The selective excitation of the v^ mode, following excitation 
of out-of-place heme doming, led to the intriguing conjecture that there may be directed 
energy transfer of the heme excitation to low frequency motions connected to backbone dis- 
placement and to protein function. The low frequency heme modes (<400 cm~^) have been 



studied using femtosecond coherence spectroscopy with a 50 fs pulse [£ 



. A series of modes 



at ~40 cm~^, ~80 cm^^, ~130 cm^^ and ~170 cm~^ were observed for several myoglobin 
derivatives. The couplings between these modes were suggested. It is a long term goal of 
our studies to understand, at the mode-specific level, how the flow of excess energy due to 
ligand dissociation leads to the selective excitation of z/4 and u-j modes. 

In our study, we ignore the transition in spin state that occurs upon ligand photodis- 
sociation and the associated electron-nuclear coupling that will no doubt be essential to 
an understanding of the "initial state" of the z/4 and z/7 vibrations following ligand pho- 
todissociation. Our focus was on the less ambitious but important question of vibrational 
energy flow on the ground state {S = 2) surface following excitation resulting from pho- 
todissociation. We employed time-dependent perturbation theory, Eq. fl50|) . to model the 
mode-specific relaxation dynamics. The initial decay process of each system mode was fitted 
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by a single-exponential function. The time constant of ~ 1.7 ps was derived for the z/4 mode 
and ^ 2.9 ps for the z/7 mode. These theoretical predictions, which make no assumptions 
regarding the VER mechanism, agree well with previous experimental results of Mizutani 



and Kitagawa for MbCO 



18| 



Vibrational energy transfer pathways were identified by calculating the 3rd order Fermi 
resonance parameters, Eq. (152|) . For the excited and v-j modes, the dominant VER 
pathways involve porphine out-of-plane motions as energy accepting doorway modes. Im- 
portantly, no direct energy transfer between the z/4 and modes was observed. Cooling of 
the five Fe-oop (Fe-out-of-plane) modes, including the functionally important heme doming 





FIG. 12: Optimized structure of FeP-Im (quintuplet S = 2 spin configuration) at the UB3LYP/6- 
31G(d) level of theory. [Reprinted with permission from Y. Zhang, H. Fujisaki, and J.E. Straub, 
J. Chem. Phys. 130, 025102 (2009). Copyright ©2009 by the American Institute of Physics.] 
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motion and Fe-Im stretching motion, takes place on the picosecond time scale. All modes 
dissipate vibrational energy through couplings, weaker or stronger, with low frequency out- 
of-plane modes involving significant imidazole ligand motion. It has been suggested that 
these couplings trigger the delocalized protein backbone motion, important for protein func- 
tion, which follows ligand dissociation in Mb. 

The 77 mode, a porphine methine wagging motion associated with Fe-oop motion, is 
believed to be directly excited following ligand photodissociation in MbCO. The coupling 
of this mode to lower frequency bath modes is predicted to be very weak. However, its 
overtone is strongly coupled to the v-j mode, forming an effective energy transfer pathway 
for relaxation on the electronic ground state and excited state surfaces. This strong coupling 
suggests a possible mechanism of excitation of the v-j mode through energy transfer from 
the 77 mode. That mechanism is distinctly different from direct excitation together with Fe- 
oop motion of the z/4 mode and supports earlier conjectures of mode-specific energy transfer 



following ligand dissociation in myoglobin 18 



V. SUMMARY AND DISCUSSION 

This chapter provides an overview of our recent work on the application of the non- 
Markovian theory of vibrational energy relaxation to a variety of systems of biomolecular 
interest, including protein backbone mimicking amide I modes in A^-methylacetamide (in 
vacuum, in water cluster, and in solvent water) amide I modes in solvated cytochrome c, 
and vibrational modes in a heme-mimicking porphyrin ligated to imidazole. We calculated 
the VER time scales and mechanisms using Eq. ( H9|) . incorporating a fiuctuating bath, 
and Eq. f l50|) . using a static bath approximation, and compared them to experiment when 
available. The theory is based on the reduced model using normal mode concepts with 3rd 
and 4th order anharmonicity, Eq. (jj]). Applying the simple time-dependent perturbation 
theory, and ensemble averaging the resulting density matrix, a non-Markovian theory of 
VER was obtained. We extended the previous theory due to Fujisaki, Zhang and Straub 



22| to more general situations where (1) the relaxing "system" has a multi-mode character 
and (2) when the system parameters depend on time (47|. We also discussed the limitations 
of our VER formulas related to the assumptions upon which the earlier theories are based. 
We are now in a position to discuss the future aspects of our work, and the connection 
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to other biolomolecular systems or nanotechnological devices. 



Relation to enzymatic reaction: The ro 



e of vibrational motions in the mechanism of 



enzymatic reactions remains controversial 93|]. In enzymology, the characterization of the 



enzymatic reaction rate is essential. Kinetic information is typically derived from substrate- 
enzyme kinetics experiments. In simulation, the free energy calculation combined with 
transition state theory is the most powerful and practical way to compute reaction rates. As 
enzymatically catalyzed reactions typically involve chemical bond breaking and formation, 
QM/MM type methods should be employed. Warshel and coworkers have examined this 
issue for several decades, and concluded that characterizing the free energy barrier is the most 
important consideration, noting that the electrostatic influence from the protein (enzyme) 
plays a key role 94] • However, Hammes-Schiffer and coworkers have identified important 






FIG. 13: Time course of the pu element of the reduced density matrix for 1/4 and exci- 
tations of f = 1. [Reprinted with permission from Y. Zhang, H. Fujisaki, and J.E. Straub, 
J. Chem. Phys. 130, 025102 (2009). Copyright ©2009 by the American Institute of Physics.] 
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situations in which VER might play a role in controlling the rate of enzymatic reactions 



93|. Furthermore, Hynes and coworkers applied the Grote-Hynes theory to the enzymatic 
reactions, and investigated the dynamic role of the environment gI- These recent studies 
indicate the importance of incorporating vibrational energy flow and dynamics as part of a 
complete understanding of enzyme kinetics. 

Relation to conformational change: The relation between vibrational excita- 
tion/relaxation and conformational change of molecules is intriguing in part because 
of the possible relation to the optimal control of molecular conformational change using 
tailored laser pulses. It is well known that there are dynamic corrections to the RRKM 
reaction rate: the simplest being 



k{E) = kRRKM{E)/{l + VR/klVR{E)) 



(53) 



where vr is the intrinsic frequency of a reaction coordinate, kivR{E) is a microcanonical 
IVR rate, and kRRKM{E) is the RRKM reaction rate 3|-l5[. Several modifications to this 
formula are summarized in [12]. It is obvious that VER affects how a molecule changes its 
shape. However, this is a "passive" role of VER. Combining RRKM theory and the local 
random matrix theory p], Leitner, Wales and coworkers theoretically studied the active 
role of vibrational excitations on conformational change of a peptide-like molecule (called 
NATMA) 95| . There are two particular modes (NH stretching) in NATMA, and they found 



that the final product depends on which vibrational mode is excited [96|. For the same 
system, Teramoto and Komatsuzaki further refined the calculation by employing ab initio 
potential energy surface 97| . A possibility to control molecular configurations of peptides or 



proteins using laser pulses should be pursued and some experimental attempts have begun 



20 



98 



99|. 



Another interesting attempt should be to address mode specific ene rgy flow associated 
with structural change. Recently Ikeguchi, Kidera, and coworkers |10Cll | developed a linear 
response theory for conformational changes of biomolecules, which is summarized in another 



chapter of this volume [37j. Though the original formulation is based on a static picture 
of the linear response theory (susceptibility), its nonequilibrium extension may be used to 
explore the relation bet ween energy flow and conformational change in proteins. In addition, 
Koyama and coworkers 10 1| devised a method based on principal component analysis for 
individual interaction energies of a peptide (and water), and found an interesting correlation 
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between the principal modes and the direction of conformational change 101 1. 

Relation to signal transduction in proteins: Though signal transduction in biology 
mainly denotes the information transfer processes carried out by a series of proteins in 
a cell, it can be interesting and useful to study the information flow in a single protein, 
which should be related to vibrational population and phase dynamics. Straub and cowork- 
ers 4l| studied such energy flow pathways in myoglibin, and found particular pathways 



rom heme to water, later confirmed exp erimentally by Ki taga wa, Mizutani and coworkers 



102| and Champion and coworkers 92|. Ota and Agard 



103l | devised a novel simulation 



protocol, anisotropic thermal diffusion, and found a partitular energy flow pathway in the 
PDZ domain protein. Importantly, the pathway they identified is located near the conserved 
amino acid region in the protein family prev iously elucidated usi ng iri formation theoretic 



approach by Lockless and Ranganathan [104f]. Sharp and Skinner [105 1 proposed an alter- 
native method, pump-probe MD, and examined the same PDZ domain protein, identifying 
alternative energy flow pat hways. Using linear response theory describing thermal diffusion, 
Ishikura and Yamato 106| discussed the energy flow pathways in photoactive yellow pro- 
tein. This method w as re cently extended to the frequency domain by Leitner and applied 



to myoglobin dimer 



1071]. Though the energy flow mentioned above occurs quite rapidly 



~ ps), there are time- resolved spectroscopic methods to detect these pathways in vitro 
20l | . Comparison between theory and experiment will help clarify the biological role of such 
energy flow in biomolecular systems. 

Exploring the role of VER in nanodevice design: Applications of the methods described 
in this chapter are not limited to biomolecular systems. As mentioned in the introduction, 
heat generation is always an issue in nanotechnology, and an understanding of VER in 
molecular devices can potentially play an important role in optimal device design. The 
estimation of thermal conductivity in such devices is a good starting point recently pursued 



by Leitner 



12| . Nit zan and coworkers studied thermal conduction in a molecular wire using 



a simplified model 108| . It will be interesting to add more molecular detail to such model 
calculations. Electr onic conduction has been one of the main topics in nanotechnology and 
mesoscopic physics 109] , and heat generation during electronic current flow is an additional 
related area of importance. 

VER in a confined environment: We have found evidence for spatially anisotropic vibra- 
tional energy flow with specific pathways determined by resonance and coupling conditions. 
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It was shown for amide I modes in cytochrome c that VER may depend on the position of 



the probing modes 



23| . making it useful for the study of inhomogeneity of th e en vironment. 



For example, an experimental study of VER in a reverse m icelle environment 110| . fullerene, 



may all be approached using 



Recently Toga shi and 



nanotube, membrane, or on atomic or molecular surfaces 
methods described in this chapter. 

Anharmonic effects in coarse-grained models of proteins: 
Mikhailov studied the conformational relaxation of elastic network models 112|. Though 
the model does not explicitly incorporate anharmonicity, small anharmonicity exists, 
resulting in interesting physical behavior relevant to biological function. Sanejouand and 
coworkers adde d exp licit anharmonicity into the elastic network models, and studied the 
ener gy storage 113| through the lens of "discrete breather" ideas from nonlinear science 



114| . Surprisingly, they found that the energy storage may occur in the active sites of 



proteins. It remains to be seen whether their conjecture will hold for all-atom models of 
the same system. 
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